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The  epicenters  of  all  89  Semipalatinsk  underground  nuclear  explosions  that  Nuttli  (1986b,  1987)  8 

studied.  57  events  for  which  the  individual  station  mb{Lg)  values  are  available  at  DARPA/NMRO 
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12/03/75)  as  well  as  about  half  of  the  Balapan  events  used  in  Nuttli's  (1986b,  1987)  study.  The 
32  missing  events  are  shown  in  diamonds. 

The  Q0  values  measured  by  Nuttli  (1986b)  with  the  coda-0  method  (top),  and  those  derived  by  16 
the  GLM  scheme  (bottom).  The  GLM  joint  inversion  gives  a  O0  of  756  for  the  Kazakh-COP 
path,  which  is  8%  higher  than  Nuttli’s  guess.  Thus  events  solely  recorded  at  the  station  COP 
(such  as  Balapan  explosion  02/10/85)  would  have  significant  discrepancy  between  Nuttli's  origi¬ 
nal  mb(Lg)  and  the  GLM/LSMF  result.  Nuttli  did  not  have  sufficient  reliable  data  for  the  path 
from  Eastern  Kazakhstan  to  COP  to  determine  the  Q0  and  C  with  the  coda-Q  method,  and 
hence  the  values  for  the  station  KON  were  used  (Nuttli,  1986b;  page  1243).  The  station  COP  is 
right  on  the  same  azimuth  as  Russian  IRIS  stations  OBN  for  which  Bennett  et  al.  (1990)  also 
find  a  O0  of  760. 
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the  attenuation  correction  in  his  mb(Lg)  determination,  therefore  the  station  residuals  (inferred 
by  LSMF)  should  be  relatively  small  for  mb(Lg) .  The  large  positive  station  term  at  COP  would 
suggest  that  Nuttli  could  have  underestimated  the  path  Q0  for  all  paths  from  Eastern  Kazakh  to 
COP.  The  bias  of  0.180  can  be  translated  to  an  increase  of  8%  in  the  path  O0,  which  would 
yield  a  Q0  estimate  identical  to  what  obtained  with  the  GLM  iterative  inversion  method. 

GLM/LSMF-reprocessed  mb(Lg )  values  versus  Nuttli’s  (1986b,  1987)  mb(Lg)  values.  The  18 

Balapan  event  02/10/85  stands  out  as  an  outlier.  For  all  other  56  events,  Nuttli’s  event  mb(Lg) 
values  appear  to  be  very  consistent  with  our  reprocessed  results. 
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SUMMARY 


This  is  our  second  semi-annual  technical  report  prepared  under  DARPA  contract 
F29601-91-C-DB23.  It  summarizes  the  research  conducted  during  April  -  October,  1992.  An 
iterative  procedure  is  developed  to  invert  for  the  path  y  and  the  event  mb{Lg)  values  simul¬ 
taneously  without  using  any  a  priori  path  y  information.  This  joint  inversion  scheme  can  be 
applied  to  data  from  multiple  test  sites  as  well.  With  this  iterative  method,  it  is  rather  easy  to 
incorporate  the  appropriate  boundary  condition  into  the  inversion.  The  iterative  scheme  is  less 
sensitive  to  rounding  errors,  and  hence  it  is  numerically  more  accurate  than  those  direct 
methods  based  on  matrix  factorization  or  Gaussian  elimination.  When  the  number  of  equa¬ 
tions  becomes  large,  the  iterative  method  is  often  the  only  practical  means  to  tackle  the  prob¬ 
lem.  The  necessity  of  incorporating  an  extra  boundary  condition  is  also  reviewed  in  this  report. 
This  joint  inversion  scheme  is  a  natural  extension  to  crustal  phases  of  the  one  we  used  in  the 
mb  analyses  ( cf .  our  first  semi-annual  technical  report  TGAL-92-05). 

The  joint  inversion  scheme  has  been  applied  to  Nuttli's  Semipalatinsk  mb(Lg)  data  set 
recovered  by  DARPA/NMRO.  The  result  indicates  that  the  majority  of  Q0  values  Nuttli  used  in 
his  pioneering  Lg  study  are  very  good  except  for  the  path  from  Eastern  Kazakh  to  COP 
(Copenhagen,  Denmark).  An  8%  increase  in  Nuttli’s  postulated  coefficient  of  anelastic 
attenuation,  Q0,  for  this  particular  path  is  necessary.  Overall,  Nuttli’s  (1986b,  1987)  mb(Lg) 
data  set  for  Semipalatinsk  explosions  appears  to  have  very  good  quality  with  a  standard  devi¬ 
ation  (of  single  observation)  of  only  0.1  magnitude  unit. 
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1 .  INTRODUCTION 


The  standard  procedure  used  in  estimating  the  source  size  of  underground  nuclear 
explosions  using  mb  measurements  has  been  to  separate  the  station  terms  from  the 
network-averaged  source  terms.  The  station  terms  thus  derived  actually  reflect  the  combina¬ 
tion  of  the  path  effect  and  the  station  effect,  when  only  those  events  in  a  close  proximity  are 
utilized.  If  worldwide  explosions  are  used  in  the  inversion,  then  the  path  effect  tends  to  be 
averaged  out  at  each  station.  In  either  case,  the  effect  due  to  the  propagation  path  alone 
would  not  be  obvious.  Recently  Jih  and  Wagner  (1991,  1992ab)  decomposed  the  station  mb 
with  the  following  joint  model: 

E(i)  +  S(j)  +  F(k(i),j)  +  e(i,j)  =  log10[A(i,j)/T(i,j)]  +  B(A(i,j))  .  [1] 

The  right-hand  side  of  [1]  is  the  conventional  raw  station  mb  of  the  i-th  event  (of  “true”  size 
E(i))  observed  at  the  j-th  station  where  A,  T,  and  B(A)  are  the  displacement  amplitude,  the 
dominant  period,  and  the  distance-correction  term,  respectively.  On  the  left-hand  side,  S 
represents  the  station  correction,  and  F(k(i),j)  is  the  path  correction  at  the  j-th  station  for 
explosions  from  the  k(i)-th  source  region.  Applying  this  joint  inversion  model  yields  a  reduction 
in  the  fluctuational  variation  of  station  mb  by  a  factor  ranging  from  1.2  to  3.0  (Jih  and 
Wagner,  1992ab). 

Underground  nuclear  explosions  at  continental  test  sites  produce  crustal  phases  that  can 
be  observed  at  near  or  regional  distances,  of  which  the  most  prominent  often  is  Lg  .  It  is 
natural  to  ask  if  the  joint  inversion  scheme  described  in  Equation  [1]  can  be  extended  to  Lg  . 
In  this  study,  an  iterative  procedure  is  presented  to  invert  for  the  path  y  and  the  event  mb(Lg) 
values  simultaneously  without  the  need  to  measure  the  path  y  with  other  methods  in  advance. 
This  joint  inversion  scheme  can  be  applied  to  data  from  multiple  test  sites  as  well.  With  this 
iterative  method,  it  is  rather  easy  to  incorporate  the  boundary  condition  into  the  inversion. 
The  iterative  scheme  is  less  sensitive  to  rounding  errors,  and  hence  it  is  numerically  more 
accurate  than  those  direct  methods  based  on  matrix  factorization  or  row  elimination.  When  the 
number  of  equations  becomes  huge,  the  iterative  method  is  often  the  only  practical  means  to 
tackle  the  problem.  The  need  of  incorporating  extra  boundary  conditions  in  this  type  of  inver¬ 
sion  problem  is  also  reviewed. 

The  iterative  joint  inversion  scheme  has  been  tested  with  Nuttli’s  Semipalatinsk  mb(Lg) 
data  set  furnished  by  DARPA/NMRO.  The  result  indicates  that  the  majority  of  Q0  values  Nuttli 
used  in  his  pioneering  Lg  study  are  very  good  except  for  the  path  from  Eastern  Kazakh  to 
COP  (Copenhagen,  Denmark).  An  8%  increase  in  Nuttli’s  postulated  coefficient  of  anelastic 
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attenuation  for  this  particular  path  is  necessary.  Overall,  Nuttli's  (1986b,  1987)  mb(Lg)  data 
for  Semipalatinsk  explosions  appear  to  have  very  good  quality  with  a  standard  deviation  (of 
single  observation)  of  0.1  magnitude  unit  [m.u.]. 


2.  JOINT  INVERSION  MODEL 

For  the  purpose  of  calculating  the  individual  mb(Lg),  Nuttli’s  (1986ab)  original  formulae 
can  be  rewritten  in  a  slightly  more  convenient  form  (Jih  and  Lynnes,  1992): 

mb{La)  ■  4.0272  +  logA(A)  +  -log(A)  +  -log[sin( - - )j  +  Y(A-1_0km) 

M  9>  a  \  ;  3  ;  2  yi  ^  1111(km/deg)^  |n(10)  ^ 

where  Y  =  Q-  ~  -T  .  Q(f)  =  Q o  ■&. 

A  is  the  epicentral  distance  in  kilometers,  A(A)  is  the  observed  0-to-peak  Lg  amplitude  meas¬ 
ured  in  the  time  domain  in  microns  [urn]  at  the  epicentral  distance  of  A  km.  The  term  —  log(A) 

3 

accounts  for  the  dispersion  when  measurements  are  made  in  the  time  domain.  The  terms 
0.5  •  log[sin(A(km)/1 11.1  (km/deg))]  and  y(A-10km)/ln(10)  correct  for  the  geometrical  spreading 
and  the  anelastic  attenuation,  respectively,  where  y  is  the  coefficient  of  anelastic  attenuation. 
For  instance,  a  seismic  source  with  1-sec  Lg  amplitude  of  110  ^im  at  10-km  epicentral  dis¬ 
tance  would  correspond  to  a  mb(Lg)  of  4.0272  +  2.0414  +  0.3333  -  1.4019  +  0.0000  =  5.000, 
same  as  what  Nuttli’s  original  formulae  would  give.  This  formula  is  appropriate  for  events  in 
eastern  North  America  and  Central  Asia. 

Now  consider  the  case  of  m  explosions  recorded  at  some  or  all  of  r.  Nations.  For  simpli¬ 
city,  we  assume  that  these  m  explosions  are  detonated  in  the  same  test  site  for  the  moment. 
The  case  of  multiple  test  sites  will  be  discussed  later.  The  linear  model  for  Lg  phases  can 
then  be  written  as 

E(i)  ~  7(j)(A(ij)-10km)/ln(10)  +  e(i.j)  =  Y(i.j)  ,  [3] 

where  Y(i,j)  =  4.0272  +  logA(i.j)  +  ~log(A(i,j))  +  ^log[sin(  — — ■  -)]  . 

3  2  1 1 1 .1  (km/deg) 

Once  the  amplitudes  and  the  locations  of  the  events  (and  hence  the  epicentral  distances,  A) 
are  available,  Y  would  be  completely  known.  Only  the  event  sizes  (E)  and  the  path-specific 
coefficients  of  anelastic  attenuation  (y)  are  the  unknown  parameters  to  be  determined.  The 
obscuring  errors  e  are  assumed  to  be  uncorrelated  and  to  uelong  to  the  same  probability  dis¬ 
tribution,  namely  a  common  Gaussian  distribution  with  zero  mean  and  variance  o2.  Our 
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model  (Equation  [3])  is  a  special  case  of  a  more  general  linear  system: 

E(i)  +  S(j)  +  e(i,j)  =  Y(i.j)  ,  for  i  =  1 . m  ;  j  =  1 . n  .  [4] 

which  can  be  expressed  in  a  matrix  formulation: 

H  X  +  e  =  Y  ,  [5] 

where  H  is  the  design  (or  observation)  matrix.  X  and  Y  are  the  column  vectors  of  unknowns 
and  observations,  respectively.  The  standard  least-squares  (LS)  solution  {viz. ,  the  one  that 
minimizes  the  residual  sum  of  squares:  RSS  a  (Y  -  HX)'  (Y  -  HX))  to  any  linear  system  with  a 
general  form  like  Equation  [5]  is 

X|_s  =  (H'H)_1H'Y  ,  [6] 

where  H'  is  the  transpose  of  H.  This  least-squares  estimator  has  many  optimality  properties. 
For  instance,  it  is  unbiased,  and  it  gives  minimum  variance  within  the  class  of  linear  unbiased 
estimators.  Furthermore,  XLS  is  also  the  Maximum-Likelihood  Estimate  (MLE)  under  the  Gaus 
sian  assumption.  It  is  straightforward  to  compute  the  uncertainty  by  using  Var[)<LS]  =  the  diag¬ 
onal  of  o^H'H)-1,  which  is  simply  scaling  the  variance  of  the  random  perturbations  by  the 
number  of  observations  associated  with  each  unknown. 


In  our  case,  however,  the  matrix  H'H  in  Equation  [6]  is  singular,  and  her  ce  the  least- 
squares  theory  can  not  be  applied  immediately  unless  the  linear  system  of  [4]  is  modified 
somewhat.  Perhaps  the  easiest  way  to  illustrate  the  indeterminacy  due  to  the  singularity  of 
the  matrix  H'H  is  that  given  any  set  of  solution  to  [4],  we  can  always  obtain  yet  another  sat  of 

solution  by  adding  a  constant  to  all  event  magnitudes,  E(i),  i  =  1 .  m,  and  subtracting  the 

same  constant  from  each  station  term  S(j)  (Jih  and  Shumway,  1989).  Alternatively,  we  can 
verify  that  the  matrix  H'H  has  zero  r  'terminant  with  linear  algebra  packages  such  as  UN¬ 
PACK,  EISPACK,  or  LAPACK.  in  this  study,  however,  a  formal  proof  is  presented  below. 
Without  loss  of  generality,  we  can  assume  that  each  of  the  m  events  is  fully  recorded  at  all  n 
stations,  then 


X  =  [E„E2,E3,  ....  Em,  S„S2 . sn]\  H'H 


n  ’^m  ^mxn 
1  nxm  nr^n 


where  lm  is  the  identity  matrix  of  order  m,  and  all  elements  of  the  m-by-n  matrix  1mxn  are  1. 
For  instance,  if  m  =  3  and  n  =  2,  then 
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H'H 


2  0  0  1  1 

0  2  0  1  1 

0  0  2  1  1 

1113  0 
1110  3 


which  is  a  doubly-bordered  band  diagonal  sparse  matrix  (Tewarson,  1973;  Press  et  a!., 
1988).  After  exactly  n  row  operations  eliminating  the  lower-left  submatrix,  1nxm,  the  deter¬ 
minant  of  H'H  can  be  computed  (up  to  a  multiplicative  constant)  as  that  of 


where  Pn  is  a  square  matrix  of  order  n  with  -^-(n-1)  on  the  diagonal  and  elsewhere.  For 
m  =  3  and  n  =  2,  (7]  becomes 


1  0  0  0.5  0.5 

0  1  0  0.5  0.5 

0  0  1  0.5  0.5 

0  0  0  1.5  -1.5 

0  0  0  -1.5  1.5 

It  suffices  to  examine  the  determinant  of  this  Pn,  or,  equivalently,  to  check  the  determinant  of 
a  square  matrix  of  order  n  with  n-1  on  the  diagonal  and  -1  elsewhere: 


n-1 

-1 

-1 

...  _i 

-1 

n-1 

-1 

...  _i 

-1 

-1 

n-1 

...  _i 

.-1 

-1 

-1 

—  n-1 

It  is  straightforward  to  prove  that,  by  mathematical  induction,  this  matrix  has  zero  determinant 
for  any  n  >  2.  Thus  the  matrix  H'H  in  our  linear  model  is  always  singular  regardless  how 
good  the  observed  amplitudes/magnitudes,  Y,  might  be.  We  therefore  need  an  extra  boun¬ 
dary  condition  to  constrain  our  linear  model  for  a  unique  solution.  The  most  commonly 
adopted  approach  in  teleseismic  magnitude  determination  is  to  have  all  station  terms  sum  up 
to  zero  (Douglas,  1966;  Blandford  and  Shumway,  1982;  Jih  and  Shumway,  1989;  Murphy 
et  al.,  1989;  Jih  and  Wagner,  1992ab;  and  many  others).  This  implies  that  larger  events 
would  tend  to  be  unchanged  whether  we  apply  the  station  corrections  or  not,  and  hence  the 
bulletin  magnitudes  of  larger  events  published  by  ISC  and  NEIC  can  be  regarded  as  more  or 
less  unbiased.  This  boundary  condition  can  be  incorporated  into  [4]  by  replacing  all  S(n)  by 

n  1 

£so).  it  not  only  reduces  the  number  of  unknowns  by  one,  but,  more  importantly,  regularizes 

i  ’ 

the  whole  linear  system  to  make  H'H  invertible.  However,  this  is  not  the  only  boundary  condi¬ 
tion  feasible.  We  can  impose  the  extra  constraint  on  E(i)  instead,  or,  even  impose  the 
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constraint  on  some  selected  stations. 


3.  ITERATIVE  INVERSION  PROCEDURE 

Once  an  extra  boundary  condition  has  been  chose",  the  inverse  matrix  of  H'H  in  Equa¬ 
tion  [6]  can  be  computed  with  matrix  factorization  (such  as  Singular  Value  Decomposition, 
[SVD])  or  Gaussian  elimination  method.  Numerical  algorithms  of  these  types  are  called  direct 
methods.  Direct  methods  can  be  impractical  if  H'H  is  large  and  sparse.  The  linear  system 
solved  in  Jih  and  Wagner  (1992ab)  for  the  optimal  mb  estimates  involved  2733  unknown 
parameters  and  24529  amplitude  measurements.  For  these  types  of  problems,  iterative 
methods  are  often  the  only  possible  method  of  solution,  as  well  as  being  faster  and  more 
accurate  than  Gaussian  elimination  and  matrix  factorization  in  many  cases.  The  largest  area 
for  the  application  of  iterative  methods  is  that  of  the  linear  systems  arising  in  the  numerical 
solution  of  partial  differential  equations.  Systems  of  orders  10,000  to  100,000  are  not  unusual 
in  aerospace  sciences,  although  the  majority  of  the  coefficients  of  the  systems  are  typically 
zeros. 

The  basic  idea  of  iterative  methods  is  that  one  starts  with  a  trial  solution  vector  X(0)  and 
carries  out  some  process  using  H,  Y,  and  X(0)  to  get  a  new  vector  X(1).  Then  one  repeats.  At 
the  k  stage,  one  uses  the  iterative  process  to  get  Xw  from  H  (or  H'H),  Y,  and  X(k-1).  The 
specific  algorithm  for  our  problem  is  summarized  in  five  steps: 


Step  0 

Set  initial  value  of  y(j)  for  j  =  1,...,  n. 

Step  1 

Compute  event  mb{Lg),  E(i),  for  i  =  1 . m: 


E(i)  =  +  Y(j)[A(i,j)-10km]/ln(10)]  , 


where  #(j)  is  the  number  of  stations  used  in  the  summation. 


Step  2 

Adjust  E(i),  i  =  1 . m,  with  the  desired  boundary  condition. 

Step  3 

Compute  the  path-specific  coefficient  of  anelastic  attenuation,  y(j),  for  j=1 . n: 

Y(j)  =  ln(10)£[  E(i)  -  Y(i,j)  ]/£[A(i,j)-10km]  . 

I  I 
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Step  4 

Fill  in  missing  paths  (i,j)  with  predicted  pseudo-observations: 

Y(i.j)  =  E(i)  +  7(j)(A(i,j)-10km)/ln(10)  . 

Step  5 

Repeat  steps  [1]-[4]  to  update  E  and  y  till  convergence. 

Although  an  arbitrary  guess  of  y  will  do  at  Step  0,  picking  an  initial  y  close  to  the  average 
across  the  whole  area  of  interest  would  speed  up  the  convergence.  The  choice  of  the  extra 
boundary  condition  (at  Step  2)  is  very  flexible.  Constraining  the  mean  of  all  event  mb(Lg ) 
seems  to  perform  extraordinarily  well,  however.  Note  that  the  pseudo-observations  predicted 
at  Step  4  are  treated  as  “good”  observations  at  steps  1  and  3  except  during  the  first  iteration 
loop. 

The  iterative  procedure  described  above  is  a  special  case  of  a  general  algorithm  known 
variously  as  the  Gaussian-Seidel  method,  the  Liebmann  process,  or  the  method  of  successive 
displacements  (Bunch  and  Rose,  1976;  Forsythe  et  al.,  1977;  Golub  and  Van  Loan,  1983; 
Spedicato,  1991).  The  major  difference  between  this  method  and  that  of  Gaussian -Jacobi 
(viz.,  the  simultaneous  displacement  method)  is  that  we  solve  for  one  component  (viz. ,  E(i)  or 
y(j))  of  the  new  vector  X  using  for  each  other  component  of  X  its  most  recently  computed 
value,  whereas  Gaussian-Jacobi  method  updates  all  unknown  parameters  simultaneously  at 
the  end  of  each  iteration  loop.  In  our  case,  Gaussian-Seidel  method  converges  faster  than 
does  Gaussian-Jacobi  method. 

The  advantage  of  our  multi-event  joint  inversion  scheme,  as  compared  to  the  simple  net¬ 
work  averaging  for  each  individual  event  that  Nuttli  (1986ab,  1987,  1988)  used,  is  that  we  can 
have  a  more  consistent  network  for  all  events.  By  "consistent"  it  means  that  the  joint  inversion 
procedure  provides  the  best  approximation  to  the  network-averaged  values  that  would  have 
been  obtained  if  all  the  events  were  recorded  at  every  station  in  the  network  (Murphy  et  al. , 
1989).  This  advantage  may  not  be  very  obvious  in  using  direct  methods  such  as  SVD  or 
Gaussian  elimination.  However,  it  would  become  quite  natural  as  revealed  by  Step  4  of  our 
iteration  scheme. 
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4.  NUTTLI’S  SEMIPALATINSK  m.  (L  )  DATA  SET 

d'  g' 

Figure  1  shows  the  epicenters  of  all  Semipalatinsk  underground  nuclear  explosions  that 
Nuttli  (1986b,  1987)  studied.  The  epicenters  used  in  this  study  (for  computing  the  epicentral 
distance,  A)  are  those  of  Lilwall  and  Farthing  (1990).  57  events  for  which  the  individual  sta¬ 
tion  mb(Lg)  values  are  available  at  DARPA/NMRO  are  listed  in  Table  1.  The  events  are 
sorted  and  identified  by  their  dates  (year,  month,  and  day)  and  the  subregion  they  are 
located.  There  are  31,  24,  and  2  events  from  Balapan,  Degelen  Mountain,  and  Murzhik  sub¬ 
sites,  respectively  (shown  in  crosses  in  Figure  1).  These  57  events  cover  most  Degelen 
events  (except  the  one  on  12/03/75)  as  well  as  about  half  of  the  Balapan  events  that  were 
used  in  Nuttli’s  (1986b,  1987)  study.  The  average  of  these  57  event  magnitudes  is  5.535. 
This  value  will  be  used  to  constrain  our  event  mb(Lg)  estimates  so  that  the  new  magnitudes 
would  be  in  alignment  with  those  of  Nuttli's  (1 986b,  1 987).  Semipalatinsk  test  range  has  long 
been  divided  into  three  test  sites  geographically,  namely,  Balapan,  Degelen,  and  Murzhik.  It 
has  also  been  suggested  by  several  recent  yield  estimation  studies  {e.g. ,  Marshall  et  al. , 
1984;  Ringdal  et  al.,  1992;  Jih  and  Wagner,  1992ab)  to  further  partition  Balapan  test  site  into 
three  smaller  subregions  for  better  results  based  on  mb  measurements.  In  this  study,  how¬ 
ever,  no  attempt  is  made  to  separate  the  explosions  according  to  their  source  media.  All 
explosions  are  assumed  to  belong  to  the  same  test  site,  and  they  share  the  same  propagation 
effect  to  any  single  recording  station. 

To  our  best  knowledge,  Nuttli’s  original  amplitude  measurements  are  no  longer  available. 
To  test  our  new  joint  inversion  scheme  (Equation  [3]),  these  157  mb(Lg)  magnitudes  are  con¬ 
verted  to  the  theoretical  ground  displacement  in  nanometers  [nm]  of  Airy  phase,  assuming  a 
constant  group  velocity  of  3.5  km/sec  and  a  period  of  1  second  using  the  following  formula: 

.  ...  «mh(l-nM-0272  ,  -0.333  ,  .  ,  .  ..  _  ,,-0.5  -1<A-10km) 

A(A)s10  9  A  ■  [sin(A/1 1 1.1)]  -exp  ,  [9] 

Nuttli's  (1986b)  path  Q0  values  ( cf .  Table  5)  are  used  in  [9]  to  determine  the  path  y,  and  the 
re-constructed  station  amplitudes  are  listed  in  Table  2.  Note  that  Nuttli  (1986b)  estimated  the 
Q0  and  C,  with  the  coda-0  method  for  10  out  of  11  WWSSN  stations,  and  he  assumed  that 
the  station  COP  (Copenhagen,  Denmark)  has  the  same  O0  of  700  as  that  of  the  station  KON 
(Konsberg,  Norway). 
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Figure  1.  The  epicenters  of  all  89  Semipalatinsk  underground  nuclear  explosions  that  Nuttli  (1986b, 
1987)  studied.  57  events  for  which  the  individual  station  mb{Lg)  values  are  available  at  DARPA/NMRO 
are  shown  in  crosses,  which  include  31 , 24,  and  2  events  from  Balapan,  Degelen  Mountain,  and  Murzhik 
subsrtes,  respectively.  These  57  events  covers  most  Degelen  events  (except  the  one  on  12/03/75)  as 
well  as  about  half  of  the  Balapan  events  used  in  Nuttli's  (1986b,  1987)  study.  The  32  missing  events 
are  shown  in  diamonds. 
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Event  * 

N 

Mean 

761123B 

3 

5.857 

7811298 

4 

6.007 

791028B 

5 

6.062 

791202B 

2 

6.050 

791223B 

6 

6.120 

800612B 

5 

5.742 

800629B 

5 

5.708 

801012B 

5 

5.918 

801214B 

4 

5.932 

801227B 

7 

5.996 

810329B 

2 

5.450 

810422B 

4 

5.968 

810913B 

6 

6.098 

811018B 

5 

6.088 

811227B 

5 

6.146 

820425B 

7 

6.126 

820704B 

2 

6.135 

821205B 

2 

6.210 

831006B 

2 

5.930 

831026B 

1 

6.100 

840219B 

1 

5.740 

840307B 

1 

5.680 

840329 B 

1 

5.970 

840425B 

2 

5.860 

840526 B 

2 

6.065 

84071 4B 

3 

6.170 

841027B 

2 

6.095 

841202B 

1 

5.970 

841216B 

1 

6.080 

841228B 

2 

6  130 

8502 10B 

1 

5980 

Table  1.  Nuttli's  Station  mb(Lg)  of  31  Balapan  Explosions  Furnished  by 


Mean  KBL  MHI  NDI  NIL  QUE  SHL  COP 


NMRO 


KON 


NUR  UME 


6.09  6.01  5.98  5.95 


6.15  _  6.10  6.16 


5.75  5.55  5.70 


5.75  5.85 


5.94 


5.96  6.06 


5.35 


5.90  5.96 


6.22 

6  14 

6.15 

6.12 

5.99  6.02 


5.94  6.01 


6.00 


6.04 


6.17  600 
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10 


11 
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5.  INVERSION  RESULTS  WITH  NUTTLi’S  DATA  SET 


We  carried  out  experiments  to  test  two  different  inversion  schemes.  In  the  first  test,  the 
re-constructed  station  Lg  amplitudes  (Table  2)  are  fed  into  the  iterative  inversion  procedure, 
with  a  boundary  condition  that  57  events  would  have  a  mean  mb{Lg)  of  5.535.  In  the  second 
experiment,  we  fed  the  157  station  mb{Lg )  values  in  Table  1  into  a  standard  LSMF  [Least 
Squares  Matrix  Factorization]  (Douglas,  1966)  inversion  package  based  on  SVD,  with  a  boun¬ 
dary  condition  that  ten  stations  (excluding  COP)  maintain  an  average  station  bias  of  zero.  For 
convenience,  these  two  experiments  are  denoted  as  "GLM”  [General  Linear  Model]  and 
“LSMF”,  respectively,  in  our  discussion. 

Note  that  Nuttli’s  station  mb(Lg)  values  have  already  been  corrected  for  the  anelastic 
attenuation.  If  all  the  Q0  values  Nuttli  (1986ab)  applied  are  appropriate,  then  his  station  Lg 
magnitudes  should  be  very  consistent.  As  a  result,  the  station  bias  should  all  be  negligibly 
small  if  a  LSMF-type  inversion  is  applied.  Conversely,  if  a  prominent  station  correction  is  still 
present,  then  this  would  mean  that  the  original  attenuation  correction  Nuttli  used  may  not  be 
quite  appropriate  for  certain  paths.  In  other  words,  our  LSMF  station  terms  will  give  a  direct 
clue  of  any  un-accounted  path  effect  such  as  that  due  to  inaccurate  attenuation  correction. 

The  GLM  inversion  scheme,  on  the  other  hand,  utilizes  the  raw  amplitude  measurements 
and  inverts  for  both  the  event  magnitudes  and  the  coefficients  of  path  attenuation  simultane¬ 
ously.  We  can  then  compare  the  resulting  path  Q0  values  with  those  obtained  with  the  coda- 
Q  method. 

Table  3  lists  the  path  attenuation  coefficients,  y  and  Q0,  inverted  with  GLM,  as  well  as 
the  corresponding  station  magnitude  corrections,  A mb{Lg),  inverted  with  LSMF.  Table  4  com¬ 
pares  the  resulting  event  mb(Lg)  values  with  those  by  Nuttli’s  (1986b,  1987).  Also  included  in 
Table  4  are  RMS  Lg  magnitudes  measured  by  Ringda!  et  al.  (1992)  using  NORSAR  and 
Grafenberg  (GRF)  recordings.  Ringdal's  Lg  magnitudes  are  noise-corrected  array  averages, 
obtained  by  applying  individual  bias  corrections  for  each  array  element. 

Figure  2  shows  the  Q0  values  measured  by  Nuttli  (1986b)  based  on  the  coda-Q  method 
and  those  derived  by  our  GLM  scheme.  Nuttli  (1986b)  initially  estimated  the  path  Q0  with  the 
coda-Q  method  of  Aki  and  Chouet  (1975)  and  Herrmann  (1980),  and  then  further  refined  his 
results  with  the  simple  network  averaging  of  the  station  magnitudes.  His  initial  and  final  Q0 
estimates  are  listed  as  O0  (1 )  and  Qo(2),  respectively,  in  Table  5.  The  GLM  joint  inversion 
gives  a  Q0  of  756  for  the  Kazakh-COP  path,  which  is  8%  higher  than  Nuttli’s  value  (Table  5). 
Thus  events  solely  recorded  at  the  station  COP  (such  as  Balapan  explosion  850210)  would 
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have  significant  discrepancy  between  Nuttli's  original  mb(Lg)  and  our  GLM  result.  Nuttli  did 
not  find  sufficient  reliable  data  for  the  path  from  Eastern  Kazakhstan  to  COP  to  determine  the 
O0  and  C  with  the  coda-0  method,  and  hence  the  value  for  the  station  KON  was  used  (Nuttli, 
1986b;  page  1243).  These  two  stations  are  about  4300-4400  km  away  from  the  Eastern 
Kazakhstan  test  site  with  7°  difference  in  the  azimuths.  The  station  COP  is  right  on  the  same 
direction  as  IRIS  station  OBN  (Obninsk,  Russia)  for  which  Bennett  et  al.  (1990)  also  find  a 
Q0  of  760. 

COP  recorded  seven  Balapan  events:  801227,  820425,  820704,  840425,  840714, 
841228,  and  850210  (viz. ,  events  marked  with  asterisk  in  Table  4).  All  these  seven  events 
have  smaller  GLM/LSMF  mb(Lg)  as  compared  to  Nuttli’s  original  network  average.  In  particu¬ 
lar,  for  Balapan  event  850210  for  which  COP  was  the  sole  source  of  Nuttli’s  Lg  recording,  the 
GLM/LSMF  mb(Lg)  is  identical  to  that  of  GRF  and  NORSAR,  while  Nuttli’s  mb(Lg)  is  0.180 
m.u.  larger.  This  discrepancy  of  0.180  m.u.  happens  to  be  the  averaged  station  bias  at  COP 
derived  with  the  LSMF  inversion  (Table  3).  Event  850210  stands  out  in  Figure  3  as  an  obvi¬ 
ous  outlier. 

Take  the  Balapan  event  841228  for  another  example,  which  Nuttli  obtained  station 
mb(Lg)  values  of  6.00  and  6.26  at  stations  NUR  and  COP,  respectively.  Subtracting  the  sta¬ 
tion  terms  of  0.004  and  0.180  would  yield  a  simple  network-averaged  mb(Lg)  of  6.038,  which 
is  closer  to  what  Ringdal  et  al.  (1992)  measured.  Balapan  event  840425  is  another  event 
recorded  at  exactly  the  same  two  stations  COP  and  NUR  with  station  mb(Lg)  of  5.85  and 
5.87,  respectively.  With  LSMF  station  corrections,  the  resulting  event  mb(Lg)  is  5.768,  which 
is,  again,  closer  to  the  GRF  mb(Lg).  There  is  no  GRF  or  NORSAR  mb(Lg)  of  820704  for 
comparison,  and  the  remaining  four  events  (801227,  820425,  840425,  and  840714)  were 
recorded  at  more  stations,  and  hence  the  biasing  effect  at  COP  tends  to  be  cancelled  out 
somewhat  by  the  network  averaging.  Nevertheless,  applying  the  station  corrections  deter¬ 
mined  with  LSMF  or  the  path  corrections  determined  with  GLM  systematically  yields  an  event 
mb(Lg)  closer  to  that  determined  with  NORSAR  and  GRF  data. 

These  11  y  values  shown  in  Table  3  have  an  average  of  0.0021  km-1.  The  5  Scandina¬ 
vian  stations  have  an  average  y  of  0.0014  km'1. 

Both  experiments  are  based  on  virtually  equivalent  linear  models  with  the  same  data  set 
as  input.  It  is  not  surprising  that  both  experiments  would  give  the  same  standard  deviation  of  e 
of  0.108  magnitude  unit,  which  is  a  measure  of  how  well  the  model  fits  the  data.  In  theory, 
this  should  imply  that  the  uncertainty  associated  with  the  mb(Lg)  of  the  same  event  should 
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remain  identical  under  these  two  methods,  even  though  the  mb(Lg)  itself  could  be  varying 
due  to  different  boundary  conditions.  However,  the  direct  inversion  result  based  on  SVD  is 
more  sensitive  the  rounding  error  of  the  computer,  and  hence  the  resulting  inverse  matrix, 
(H'H)-\  may  not  be  as  accurate  as  that  obtained  with  the  iterative  inversion  method.  This  is 
indicated  by  the  slight  deviation  from  the  theoretically  predicted  uncertainty  of  o/V#(j)  in  some 
of  the  LSMF-esti mated  standard  errors  in  Table  4.  Note  that,  however,  both  methods  give 
almost  identical  event  mb(Lg)  values,  even  though  two  seemingly  irrelevant  boundary  condi¬ 
tions  were  used.  If,  a  more  popular  boundary  condition  that  the  1 1  station  terms  have  zero 
mean  were  used  instead,  then  each  station  term  {viz. ,  A mb(Lg)  in  Table  3  or  S(j)  in  Equation 
[4])  would  be  reduced  by  0.016  m.u.,  and  each  LSMF  mb{Lg)  in  Table  4  would  then  be 
increased  by  0.016  m.u.  accordingly.  This  will  not  affect  our  observation  that  the  station  resi¬ 
dual  at  COP  (now  0.164  m.u.)  appears  to  be  too  large. 


Table  3.  GLM-inverted  Coefficient  of  Anelastic  Attenuation  versus  LSMF-inferred  Station  Corrections 

Station 

GLM  Attenuation  Coefficient 

LSMF  S(j) 

Code 

m 

N° 

Paths 

_ 

y  (1/km) 

Qo 

Am^Lg) 

■1  ■■ 

WM  H 

COP 

12.433 

55.683 

7 

0.00119 

756 

+0  180±0.051 

KBL 

69.043 

34.541 

15 

000253 

355 

-0.028+0.040 

KEV 

27.007 

69.755 

1 1 

0.00162 

553 

-0  116±0.037 

KON 

09.598 

59  649 

9 

0  00126 

713 

+0  045±0  039 

MHI 

59.494 

36.300 

10 

0  00240 

374 

-0.038 ±0  045 

NDI 

77.217 

28683 

18 

0.00295 

305 

+0.049±0  028 

NIL 

73.252 

33  650 

21 

0  00244 

367 

-0  065±0  034 

NUR 

24.651 

60  509 

34 

0  00154 

581 

+0  004+0  023 

QUE 

66.950 

30  188 

7 

0  00296 

303 

+0  033±0  042 

SHL 

91.883 

25  567 

15 

0  00257 

349 

+0  089±0  030 

UME 

20  237 

63  815 

10 

_ 

0  00143 

627 

+  0  026+0  036 

15 


Figure  2.  The  Q0  values  measured  by  Nuttli  (1986b)  with  the  coda -O  method  (top),  and  those  derived 
by  the  GIM  scheme  (bottom).  The  GLM  joint  inversion  gives  a  O0  ol  756  lor  the  Kazakh-COP  path, 
which  is  8%  higher  than  Nuttli's  guess.  Thus  events  solely  recorded  at  the  station  COP  (such  as 
Balapan  explosion  02/10/85)  would  have  significant  discrepancy  between  Nuttli's  original  mb(Lg)  and  the 
GLM/LSMF  result.  Nuttli  did  not  have  sufficient  reliable  data  for  the  path  from  Eastern  Kazakhstan  to 
COP  to  determine  the  Q0  and  (,  with  the  coda-0  method,  and  hence  the  values  for  the  station  KON 
were  used  (Nuttli,  1986b;  page  1243)  The  station  COP  is  right  on  the  same  azimuth  as  Russian  IRIS 
stations  OBN  for  which  Bennett  el  al  (1990)  also  find  a  O0  of  760 


16 


0.180 

X  0.090 

X 

0.060 

X 

0.045 

X 

0.036 

X 

0.030 

X 

0.026 

X 

0.023 

X 

0.020 

X 

0.018 

o 

-0.018 

o 

-0.020 

o 

0.023 

o 

-0.026 

o 

-0.030 

o 

-0.036 

o 

-0.045 

0 

-0.060 

o 

-0.090 

T) 

-0  ISO 

Figure  3.  The  mb(Lg)  station  terms  derived  with  LSMF  inversion.  Nuttli  (1986b,  1987)  already  included 
the  attenuation  correction  in  his  mb{Lg)  determination,  therefore  the  station  residuals  (inferred  by  LSMF) 
should  be  relatively  small  for  mb(Lg) .  The  large  positive  station  term  at  COP  would  suggest  that  Nuttli 
could  have  underestimated  the  path  O0  for  all  paths  from  Eastern  Kazakh  to  COP.  The  bias  of  0.180 
can  be  translated  to  an  increase  of  8%  in  the  path  O0,  which  would  yield  a  Q0  estimate  identical  to 
what  obtained  with  the  GLM  iterative  inversion  method 


Lg[GLM 1  vs.  LgfNuttli ] 


mb(Lg[GLM]) 


Figure  4.  GLM/LSMF-reprocessed  mb{Lg)  values  versus  Nuttlis  (1986b,  1987)  mb(Lg)  values.  The 
Balapan  event  02/10/85  stands  out  as  an  outlier.  For  all  other  56  events,  Nuttli’s  event  mb{Lg)  values 
appear  to  be  very  consistent  with  our  reprocessed  results. 


18 


Table  4.  Comparison  of  Various  mb(Lg) 


Event 

RMS  Lg  (Ringdal  et  al. ,  1992) 

Nuttli's  mb(Lg) 

GLM  mb(Lg) 

LSMF  mb(Lg) 

Date 

NORSAR 

GRF 

WWSSN 

WWSSN 

WWSSN 

761 1238 

5.792 

5.85710.055  3 

5  82310.062 

5820+0  065 

781129B 

5.973 

5.886 

6.00710.030  4 

5  98210.054 

5  981+0  056 

791028B 

6.053 

6.046 

6.062±0.013  5 

6.02110.048 

6.01910.050 

791202B 

5.917 

5.938 

6.05010.030  2 

602610076 

602510  080 

791223B 

6.045 

6,12010.017  6 

6  08010.044 

6079+0.046 

80061 2B 

_ 

5.571 

5.74210.059  5 

5.745+0.048 

5.743+0.050 

800629B 

5.683 

5.746 

5.70810.074  5 

5.714+0.048 

5.714+0.050 

801012B 

5.925 

5.933 

5.91810.036  5 

5.915+0.048 

5.91510.050 

801214B 

5.929 

5.944 

5.93210.035  4 

5.93110.054 

5.93210.056 

801227B * 

5.939 

5.885 

5.996+0.038  7 

5.95910.041 

5.95810.042 

810329B 

5.556 

5.437 

5.45010.100  2 

5.498+0.076 

5.497+0.083 

810422B 

5.908 

5.956 

5  96810.025  4 

5.96910  054 

5.968+0.056 

810913B 

6.113 

6.106 

6  09810.065  6 

6  08310.044 

6.08210.046 

811018B 

5.985 

5.956 

6.08810.066  5 

6  07510  048 

6.07410.050 

811227B 

6.074 

6.092 

6.14610  067  5 

6.13310.048 

6.13210.050 

820425B  * 

6.077 

6.058 

6  126+0.056  7 

6.087+0.041 

6.086+0.042 

820704B  * 

6.13510.015  2 

6.02310.076 

6022+0083 

821205B 

5.990 

6.002 

6.21010.090  2 

6  165+0.076 

6.16310.079 

831006B 

5.870 

5.843 

5  93010070  2 

5.98610076 

5.986+0.080 

831026B 

6.000 

6.036 

6.1001  1 

6.09610  108 

6.09610.111 

8402 19B 

5.725 

5.7401  1 

5.735+0  108 

5.736+0.111 

840307B 

5.698 

5.6801  1 

5.67710.108 

5.676+0.111 

840329B 

5897 

5.957 

5.9701  1 

5  965+0.108 

5.966+0.111 

840425B * 

5.870 

5.803 

5  86010.010  2 

5  76810  076 

5  76810.081 

840526B 

6.072 

6  128 

6.06510.045  2 

6.12110  076 

6  12110  080 

84071 4B  * 

6.054 

6.064 

6.17010.208  3 

6.148+0.062 

6  14710  066 

841027B 

6.085 

6.145 

6.09510.055  2 

6  151+0  076 

6.15110.080 

841202B 

5.880 

5.860 

5.9701  1 

5.966+0.108 

5.96610.111 

841216B 

6.048 

6.038 

6  0801  1 

6.07610  108 

6.076+0  111 

841228B  • 

5  985 

5.947 

6  130+0  130  2 

6  038+0  076 

6  038+0  081 

850210B  * 

5  803 

5  801 

5980+  1 

5  801+0  108 

5  80010  120 

events  recorded  at  COP 
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RMS  Lg  (Ringdai 


NORSAR 


Table  4  Comparison  of  Various  mB(Lg) 


et  at. ,  1992)  Nuttli's  rnb(Lg) 


WWSSN 


GLM  mb(Lg) 

LSMF  mb(Lg) 

WWSSN 

WWSSN 

5.118+0.076 

5  1 1 710.083 

4.975+0.108 

4.97510.113 

5.172+0.076 

5.17210.083 

4.90710.076 

4  90710.083 

4  82610.076 

4.82710.083 

4.87610.108 

4.87510.113 

4.57910.076 

4  57710.083 

4.47710.076 

4.47710.083 

5.42410.062 

5.42510.066 

5.36610.076 

5.36710.082 

5.06210.108 

5.05110.112 

5.167+0.076 

5.16710.083 

4.98610.076 

4  98710  083 

4.92910.076 

4  93310.079 

4.66910.108 

4.66810.115 

4.83410.108 

4.83810.117 

5.33810.054 

5.338+0.056 

4.964+0.076 

4.96910.081 

4.78510.062 

4  784+0.065 

4.62010.108 

4.62810  117 

5.429+0.048 

5.428+0.050 

5.045+0.054 

5.054+0.056 

5  33310.062 

5333+0.065 

5.698+0.054 

5.698+0.056 

5  11110  076 

5  111+0.079 

505010076 

5.052+0.082 

7307100 


731026D 


7401 30b  D 


740516D 


74071 00 


740913D 


741216aD 


741216bD 


750220D 


7503110 


7506080 


7508070 


7601 15D 


760421aD 


760723D 


761230D 


770329D 


770425D 


770730D 


780326D 


780422D 


780728D 


7811290 


730419M 


780319M 


5.07010.070  2 


4.9101 _ 1 


5.125±0.075  2 


4.860±0.000  2 


4.780+0.090  2 


4.810+ _ 1 


4.530+0,020  2 


4.430+0.050  2 


5.410+0.058  3 


5.31510.055  2 


5.100± _ 1 


5.120+0.120  2 


4  940+0.010  2 


4.96010.020  2 


4.6401 _ 1 


4  8001 _ 1 


5.33510.018  4 


4.975±0.085  2 


4.763+0.054  3 


4.5901 _ 1 


5.436+0.069  5 


5.08010.089  4 


5.33010.023  3 


5.740+0.046  4 


5.080+0.020  2 


5.000+0.050  2 


:  from  Ringdai  (personal  communication.  1991). 


We  can  translate  each  LSMF  station  residual  in  Table  3  to  a  new  estimate  of  the 
coefficient  of  anelastic  attenuation  with  the  following  formula: 

Q0(new)  =  — - — - - - 

1  _  *mb(Lg)  •  ln(  1 0)  •  U  [8] 

Oo(°lcl)  (A-IOkm)  •  71 

This  is  essentially  the  same  approach  Nuttli  (1986ab,  1987,  1988)  and  Patton  (1988)  used 
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previously.  The  resulting  new  Q0  values  are  listed  in  Table  5  as  Q0(4).  Excellent  agreement 
is  obtained  when  comparing  Q0(4)  t0  those  Q0{3)  derived  with  GLM.  Once  again,  this 
confirms  that  Q0  for  Kazakh-COP  path  should  be  higher  than  Nuttli’s  guess  of  700.  For  other 
paths,  however,  Nuttli’s  original  estimates  based  on  the  coda-0  method  are  generally  rather 
consistent  with  our  results  based  on  the  Lg  -amplitude  attenuation. 


Table  5.  Revised  Qq  for  Semipalatinsk-WWSSN  Paths 

Station 

• 

A 

Nuttli  (1986b) 

This  Study 

Code 

(km) 

c 

Qo(i' 

Qo  (2) 

Qy  3) 

Q0(4).  A  mb{Lg) 

COP 

4350 

(0.4) 

— 

(700) 

756 

756,  +0  180 

KBL 

1900 

D 

360 

336 

355 

332.  -0.028 

KEV 

3475 

0.4 

580 

554 

553 

529,  -0.116 

KON 

4375 

0.4 

700 

700 

713 

713,  +0.045 

MHI 

2150 

0.5 

380 

380 

374 

374,  -0.038 

NDI 

2375 

0.6 

300 

312 

305 

317,  +0  049 

NIL 

1875 

0.6 

380 

354 

367 

343,  -0.065 

NUR 

3525 

0.4 

580 

580 

581 

581,  +0.004 

QUE 

2425 

0.6 

300 

318 

303 

322,  +0.033 

SHL 

2925 

B 

300 

340 

349 

349,  +0.089 

UME 

3750 

0.4 

620 

591 

627 

597,  +0.026 

rounded  to  the  nearest  25  km  (from  Table  1  of  Nuttli.  1986b) 

Q0:  (1)  Nuttli's  initial  measurement  based  on  the  coda-Q  method;  (2)  Nuttli's  refined  Q0  based  on  the  single  event  network  averaging.  (3)  GLM  result;  (4)  LSMF 
result 
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6.  DISCUSSION  AND  CONCLUSIONS 

The  GLM  algorithm  presented  in  this  study  can  determine  the  path  y  as  well  as  the  event 
mb(Lg)  values  simultaneously  without  the  need  to  measure  the  path  y  with  other  methods  in 
advance.  A  staged  LSMF  inversion  scheme  has  also  been  tested  in  this  study.  Both  GLM 
and  LSMF  experiments  yield  very  consistent  results  suggesting  that  the  majority  of  O0  values 
Nuttli  used  in  his  pioneering  Lg  study  of  Semipalatinsk  explosions  (Nuttli,  1986b,  1987)  are 
very  good  except  for  the  path  from  Semipalatinsk  to  COP.  Nuttli’s  biased  attenuation  correc¬ 
tion  at  COP  resulted  in  overestimated  mb{Lg)  for  7  Balapan  events  relative  to  the 
GLM/LSMF-recomputed  event  mb(Lg)  values.  In  the  extreme  case,  his  mb(Lg)  estimate  for 
the  Balapan  event  of  10  Feb  85  (which  was  determined  solely  based  on  the  COP  recording 
alone)  is  biased  high  by  about  0.18  m.u.  An  8%  increase  in  the  postulated  coefficient  of  ane- 
lastic  attenuation  for  this  particular  path  is  necessary.  Other  than  this  problem,  Nuttli’s  (1986b, 
1987)  mb(Lg)  data  set  for  Semipalatinsk  explosions  has  very  good  quality  with  a  standard 
deviation  (of  single  observation)  of  0.1  m.u. 

Our  GLM  algorithm  is  implemented  with  an  iterative  procedure  instead  of  utilizing  the 
standard  matrix  factorization  or  row  elimination.  The  advantages  of  the  iterative  inversion 
scheme  are  obvious: 

[1]  It  is  rather  easy  to  incorporate  the  boundary  condition  into  the  inversion  scheme  without 
altering  the  simple  equations  describing  the  underlying  physical  model.  That  is,  the  matrix  ele¬ 
ments  of  H  will  remain  unchanged. 

[2]  The  iterative  scheme  is  numerically  more  accurate  simply  because  it  is  less  sensitive  to 
rounding  errors. 

[3]  When  the  number  of  equations  becomes  huge,  the  iterative  method  is  often  the  only  prac¬ 
tical  means  to  tackle  the  problem. 

The  only  possible  disadvantage  is  that  the  implementation  of  iterative  methods  usually 
depends  on  the  actual  application  field.  Therefore  it  may  not  be  always  easy  to  find  a 
general-purpose  linear-algebra  library  ready  for  use.  On  the  other  hand,  packages  suitable  for 
direct  methods  are  available  in  almost  every  commercial  software  library.  However,  this 
should  not  be  considered  as  a  serious  drawback  or  limitation  of  the  iterative  method.  Imple¬ 
menting  the  iterative  procedure  presented  in  this  study  is  just  as  simple  as  coding  up  a  pro¬ 
gram  to  call  the  standard  linear  algebra  library. 

The  boundary  condition  used  in  our  GLM  experiment  was  that  the  resulting  event 
mb(Lg )  values  would  have  the  same  mean  as  that  of  Nuttli’s  original  mb{Lg)  values.  Magni- 
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tudes  based  on  the  short-period  teleseismic  recordings  corrected  for  the  station  amplifications 
as  well  as  the  path  effects  (due  to  the  focusing/defocusing  near  the  source  region  etc.)  such 
as  those  in  Jih  and  Wagner  (1992ab)  may  provide  a  good  constraint  as  well.  Many  Lg  studies 
have  used  the  ISC  mb  for  “normalizing"  their  Lg  magnitude  scale  (e.g. ,  page  2146  of  Nuttli, 
1986a;  page  128  of  Israelson,  1992). 

If  the  test  data  all  originated  from  the  same  test  site,  it  would  be  very  difficult  to  further 
decompose  the  path  corrections  into  the  station  amplification  (due  to  the  local  geology)  and 
the  propagation  path  effect.  If,  however,  data  from  multiple  test  sites  become  available,  then 
it  is  rather  straightforward  to  adapt  our  GLM  inversion  model  to  a  more  general  formulation 
which  is  almost  identical  to  Equation  [1]  (as  used  in  Jih  and  Wagner  (1992ab)): 

E(i)  +  SO)  -  'V(k(i),j)(A(i,j) — 1 0km)/ln(1 0)  +  e(i,j)  =  Y(i,j) .  [3*] 

where  k(i)  represents  the  k-th  test  site  where  the  i-th  explosion  is  detonated. 

Israelson  (1992)  applied  the  standard  LSMF  method  to  RMS  Lg  magnitudes  measured 
on  hand-digitized  Soviet  analog  waveforms,  and  he  obtained  a  suite  of  large  station  correc¬ 
tions  ranging  from  -0.309  m.u.  at  NRI  (Norilsk,  Siberia)  to  0.664  m.u.  at  UZH  (Uzhgorod),  on 
top  of  an  assumed  average  y  of  0.0012km-1  for  Russian  Platform.  Israelson  (1992)  finds  that 
his  station  corrections  seem  to  increase  with  the  local  ambient  noise  level  at  the  receiver 
sites.  There  are  other  mechanisms  that  could  reduce  or  enhance  the  Lg  amplitude  level,  and 
yet  they  ase  not  relevant  to  the  local  amplification.  The  blockage  of  Novaya  Zemlya  Lg  as 
observed  at  KEV  due  to  the  thick  sedimentary  layer  in  Barents  Sea  is  a  good  example 
(Baumgardt,  1991).  Whatever  mechanism  it  might  actually  be,  for  explosions  confined  in  a 
small  source  region,  either  applying  the  extra  bias  corrections  at  the  receivers  (as  did  in  Isra¬ 
elson,  1992;  and  in  many  mb  and  Ms  studies)  or  applying  the  path-wise  y  correction  (as  did 
in  the  work  of  Nuttli,  Patton,  and  many  others)  would  serve  equally  well  for  the  purpose  of 
event  magnitude  determination.  However,  if  the  seismic  events  spread  over  a  vast  area  so 
that  they  can  no  longer  be  grouped  into  just  a  few  isolated  test  sites,  then  maintaining  a  tabu¬ 
lar  source-station  path  corrections  may  not  be  the  most  convenient  approach.  In  such  a  case, 
a  2-dimensional  map  to  describe  the  Lg  Q  variation  would  be  necessary.  The  individual  path 
y  values  determined  with  GLM  or  the  coda-0  method  can  and  should  be  utilized  to  establish 
the  2-dimensional  (i.e.,  regionalized)  Q  map  with  the  tomographic  inversion.  Establishing  the 
regionalized  O  map  in  advance  with  a  suite  of  carefully  measured  events  for  seismotectonic 
regions  of  proliferation  concern  is  particularly  important,  since  there  may  not  be  sufficient  data 
from  the  suspected  clandestine  tests  for  the  coda  O  determination  (Jih  and  Lynnes,  1992). 
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